Measuring kinetic coefficients by molecular dynamics simulation of zone melting 
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Molecular dynamics simulations are performed to measure the kinetic coefficient at the solid-liquid 
interface in pure gold. Results are obtained for the (111), (100) and (110) orientations. Both Au(100) 
and Au(110) are in reasonable agreement with the law proposed for collision-limited growth. For 
Au(lll), stacking fault domains form, as first reported by Burke, Broughton and Gilmer [J. Chem. 
Phys. 89, 1030 (1988)]. The consequence on the kinetics of this interface is dramatic: the measured 
kinetic coefficient is three times smaller than that predicted by collision-limited growth. Finally, 
crystallization and melting are found to be always asymmetrical but here again the effect is much 
more pronounced for the (111) orientation. 
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I. INTRODUCTION 

Solidification of pure elements is of technological in- 
terest because the way a given material solidifies usu- 
ally affects its structure and, as a consequence, its final 
elastic and other macroscopic properties. From a fun- 
damental point of view, interest in free and directed so- 
lidication comes from the underlying nonlinear physics, 
morphological instabilities being at the origin of generic 
microstructures such as dendrites or cells. 

Important theoretical and numerical contributions 
have been made to solve this difficult physical problem 
jjj . Recently, a quantitative phase field model was intro- 
duced A subsequent refinement, consisting in solving 
the diffusion equation with the help of Brownian walkers, 
permitted to bridge the wide gap between the capillary 
and diffusion lengths, allowing direct comparison with 
experiments ||. As a consequence, there is currently 
an increasing need for accurate values of the interfac re- 
sponse functions that are used as input parameters for 
realistic phase field simulations. 

In the case of a pure element, the surface tension j£ mn 
must be known as a function of the interface orientation 
(£mn). In addition, the kinetic coefficient A%mi(Ii) giv- 
ing the relation between the interface velocity and the 
interface temperature Ti, should also be known for the 
different orientations. For a binary alloy, temperature de- 
pendence of the solute diffusion coefficient, D(T), as well 
as velocity and orientation dependence of the segregation 
coefficient k£ mn (Vi) are also necessary. 

Both k and \x are hardly accessible in the experiments 
and convection effects often lead to overestimated val- 
ues of diffusion coefficients. Different simulation schemes 
have thus been proposed as an alternative. Such nu- 
merical experiments have been rendered possible by the 
discovery of realistic interatomic potential models, such 
as, in the case of metals, the embedded atom model 



(EAM) D, the glue model (GM) f§ and the effective 
medium theory (EMT) B. In the near future, the in- 
crease of computer power should open the possibility to 
address the case of more complicated materials like semi- 
conductors, molecular crystals and organic compounds, 
for which potentials do not simply reduce to pair inter- 
actions. Very recently, the functions 7£ mrl and /J,e mn (Ti) 
have been determined and used in phase field simula- 
tions of dendritic growth for pure nickel [Q. The good 
quantitative agreement found between experiments and 
simulations is promising and should stimulate in the near 
future the construction of other material-dedicated phase 
field models. 

New methods for the determination of the func- 
tions jimn and ki mn (Vi) have been recently proposed 
In the present paper, we rather concentrate on 
IJ-£mn(Ti). The kinetic response of a solid- liquid inter- 
face has been simulated quantitatively in the 80's by 
Broughton, Gilmer and Jackson (BGJ) for a Lennard- 
Jones (LJ) potential and a (100) orientation |10|. These 
authors showed that growth is not diffusion-limited but 
rather that the interface velocity is related to the mean 
kinetic energy of the atoms. For this collision-limited 
growth regime, the growth rate should be directely pro- 
portional to the distance between two successive layers 
dtmn- Indeed, since the liquid atoms do not diffuse to 
choose their adsorption sites but are almost instanta- 
neously incorporated into the solid, the larger di mn the 
more effective and the faster the advance of the solid- 
liquid interface should be. The analytical expression for 
the growth velocity of a rough solid-liquid interface reads 
ill: 
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dimn being the interplane spacing, lS.fl the chemical po- 
tential difference between solid and liquid phases, T the 
absolute temperature, k the Boltzmann constant, and V 
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the thermal velocity. This law is confirmed by molecu- 
lar dynamics simulations for the (100) and (110) orien- 
tations: the expected y2 ratio between the correspond- 
ing kinetic coefficients is well recovered for several metals 
cristalizing in a face centered cubic (fee) structure (Ni, 
Ag and Au) |p2 13 1 . Nevertheless, for these rough ma- 



terials, growth of the (111) interface does not obey this 
simple law: according to Eq. (1), the (111) orientation 
should be much faster, and what is found is precisely the 
opposite. Burke, Broughton, and Gilmer |ll]] attribute 
this slowing-down to the growth of competing fee and 
hep domains in the solidfying layer, followed by the elim- 
ination of the defect lines between the two phases. 

Another question associated with solid-liquid inter- 
faces is that of symmetry between solidification and melt- 
ing kinetics. Asymmetry has been already observed in 
different systems. It is not really surprising for faceted 
materials like silicium where solidification involves nuclc- 
ation while melting does not. The question is more del- 
icate when one considers rough materials with collision- 
limited growth. Indeed, available results are controver- 
sial: if asymmetry has been found for a Na(100) interface 
Q, it has not been observed for a LJ(100) More 
surprisingly, in the latter case an opposite asymmetry 
(crystal growing faster than the melt) can be found, de- 
pending on the way the solid germ is prepared. 

In this paper, we address the above questions concern- 
ing the growth of a rough solid-liquid interface. We first 
present our implementation of a non-equilibrium molecu- 
lar dynamics scheme for a zone melting experiment. The 
second section is devoted to the study of (100) and (110) 
orientations. The special case of (111) growth is exam- 
ined in section III and asymmetry between melting and 
solidification in section IV. Finally, a summary of the 
different results and a discussion are given in the last 
section. 



II. SIMULATION PROCEDURE 

For this study we use the Ercolessi glue potential for 
Au Q . In this formalism the total potential energy for a 
system of N atoms is given by: 
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The first term is a classical pair interaction. In the second 
term, rii is the coordination of atom i, 



N 
3 = 1 
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where p(rij) is a function of the interatomic distance ry, 
with a cut-off radius of 3.9A here. The energy function 



FIG. 1. A typical simulation box with periodic bound- 
ary conditions in all three directions ((HI) solid-liquid inter- 
faces). Atoms in dark grey are within the hot and cold slices 
where temperature is fixed. 



U is the glue term associating an extra potential energy 
to atom i as a function of its coordination. This glue 
potential has demonstated its efficiency in predicting the 
physical properties of gold as well as in describing sev- 
eral experimentally observed phenomena such as surface 
melting and surface reconstructions jlfj. 

A distinctive feature of our method is to simulate a 
zone melting experiment in which both a solidification 
and a melting front are simultaneously advancing at a 
fixed velocity V. This velocity is that of the virtual fur- 
nace which imposes two symmetric thermal gradients. 
The particle coordinates are defined in a reference frame 
moving at velocity V in the z direction, so that after 
equilibration the positions of the two interfaces are fixed 
in the simulation box. Heat transport from the furnace 
is simulated by imposing one temperature below and one 
above the melting point inside two distant slices, 20A 
each in thickness (Fig. 1). Within each slice, tempera- 
ture is kept constant by using a classical velocity rescaling 
procedure Q . Periodic boundary conditions are applied 
in the three directions. More details about the numerical 
method can be found in a recent study of solute trapping 
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FIG. 2. Temperature and energy profiles along the z axis 
perpendicular to the interfaces. 
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FIG. 3. Caloric curves for V = (circles) and V — 15 
ms _1 (diamonds). For non zero velocity, the kinetic effects 
split the curve in two parts. The dotted and full straight 
lines represent respectively the functions Ei(T), El(T) and 
Es(T). 



in a LJ binary alloy, where a similar simulation technique 
was used ||. 

First, the fee solid and the liquid are equilibrated 
separatly at zero pressure and at a temperature close 
of the melting point. Our smallest system has a size 
S ^ 20 x 20A 2 in cross-section, that is about 64 atoms 
per layer. After equilibration, the solid and liquid are 
brought into contact and plunged in the temperature gra- 
dient imposed by the two temperature-controled slices. 
The total system is about 220 A in height. After a sec- 
ond equilibration period (during which the velocity of 
the furnace is zero), the two interfaces reach a station- 
ary position and we roughly have 50% of solid and liquid 
(see Fig. 1). Fig. 2 shows the temperature and energy 
profiles along the z axis. 

Combining the two profiles to eliminate the z coordi- 
nate, one obtains a caloric curve, i.e., a plot of energy as 
a function of temperature. In Fig. 3, the caloric curves 
obtained for two different values of the pulling velocity V 
are displayed. For V = 0, the data points corresponding 
to the solidification and the melting fronts merge onto 
the same curve: no kinetic effects are at play and the 
interface temperature is the equilibrium melting temper- 
ature To ~ 1330K. When a velocity is imposed, a dy- 
namical hysteresis appears on the caloric curve. Kinetic 
effects split the curve in two distinct parts: the interface 
temperature of the solidification front decreases while it 
increases on the melting front. We can deduce both in- 
terface undercooling and interface superheating from this 
plot. An interest of this method is that, as said before, 
the interface is fixed in the reference frame of the simu- 
lation box, so that statistics are easy to record. A tipical 



run lasts 10 6 MD steps (3.5 x 10 3 ps), so that the atoms 
in the system solidify and melt several times. According 
to a recent study by Tepper et al. |l5[ ], we know that 
the melting kinetics can be affected by the way the solid 
is equilibrated. Thus multi-cycling is necessary to mimic 
the melting of a real solid, usually resulting from previous 
solidification(s). 

To conclude this section, the method used to estimate 
the interface temperature T from the caloric curves is 
described. We assume the energy of atoms lying at the 
interace, E^, to be a weighted average of the perfect solid 
and liquid energies at the same temperature T. 

E i {T)=aE s (T i ) + (l-a)E L (T i ). (4) 

Linear relations, E S (T) = a s T + b s , and E L {T) = a L T + 
bh-, are fitted to the data points obtained on the low and 
high temperature side, respectively (Fig. 3). The curve 
Ei(T) is thus a line with a slope 

p = aa s + (1 - a)a L . (5) 

The value of coefficient a is then extracted from the 
caloric curve at zero velocity, for which Tj must be equal 
to To (Fig. 3). Finally, the interface temperature is given 
by the intersection of the line Ei(T) with the caloric 
curve. An alternative method consists in building an or- 
der parameter that distinguishes between solid and liquid 
atoms [§,[l8) : a plot of this order parameter as a function 
of temperature also gives an interface temperature. We 
have checked that the two methods give equivalent re- 
sults. 
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III. GROWTH OF (100) AND (110) INTERFACES 

In this section we compute the kinetic coeficient for 
the Au(100) and Au(llO) interfaces, using the method 
described above. We concentrate here on pulling veloc- 
ities ranging between V = 5 ms" 1 and V = 30 ms™ 1 , 
for which kinetics remain linear. We also perform a 
few simulations at higher velocities, where kinetics de- 
viates from linearity, but comments on nonlinear effects 
are postponed to the concluding section. In Fig. 4, we 
plot the interface velocity as a function of the measured 
undercooling To — T,-. Linear fits to the law 

V = Ht mn (To — Ti) (6) 

give the following estimates for the two kinetic coeffi- 
cients: 

(4 00 =23.1 ±1.0 cms^K" 1 (7) 

and 

l4 w = 15.5 ±1.0 cms^K- 1 . (8) 

However, finite-size effects are expected to bias these 
estimates because the system cross-section area, So 
20 x 20 A 2 , is rather small. 

Additionnal runs are thus performed in order to quan- 
titatively evaluate finite-size effects. The pulling velocity 
is fixed to V = 15ms -1 , the system height to H ~ 222A, 
and the cross-section area S is progressively increased. 
We define the normalized kinetic coefficient Hn(S) as the 
ratio of the kinetic coefficient obtained at size S to that 
obtained at size S = So [Eqs (7-8)]. As shown in Fig. 5, 
the size effects are important and the kinetic coefficients 
appear to converge only for S ~ 100 x 100A 2 . For the 
(100) direction, there is a decrease of about 20 percent 



and we obtain roughly the same behavior for the (110) 
interface. This size effect has never been reported in the 
past for (100) and (110) orientations: the fact that Hoyt 
and co-workers do no find size effects for these two orien- 
tations || is certainly due to the fact that their smaller 
system is larger than ours. We can now propose extrap- 
olated values for the kinetic coefficients: 

/iioo = 18.8 ± 1.0 cms^KT 1 (9) 

Huo = 12.6 ± 1.0 cms^KT 1 (10) 

The corresponding ratio ^ioo/mho — 1-49 ± 0.15 is in 
good agreement with the value y/2 predicted by Eq. 
(1). Hence, the assumption of collision-limited growth 
for (100) and (110) orientations is confirmed to be the 
relevant one. At this point, we can compare our results 
with those of Hoyt et al. for gold [H| . If they also find 
a V2 ratio between their two orientations, their /i val- 
ues are larger than ours by a factor 1.8. Linearizing the 
expression given by BGJ, we find 

V-T^Tr^pa-Ti) (11) 

for the interface velocity. The potential used by Hoyt 
et al. gives a melting point Tq of 1090K [^9| much 
smaller than the value 1330K obtained with Eroclessi po- 
tential. Introducing this temperature shift in Eq. (11) 
roughly accounts for the discrepancy between the values 
of fi. Since Ercolcssi potential gives a melting point much 
closer to the experimental one, it should be also the case 
for our estimates of the kinetic coefficients. 

In order to understand the origin of the size effects on 
the value of the kinetic coefficient, we take now a closer 
look at the in-plane structure of gold layers in the vicin- 
ity of the solid-liquid interface. We compute a density 
profile along the z axis from which we are able to sepa- 
rate atoms belonging to different layers. Deep in the solid 
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the in-plane square structure of the (100) orientation is 
effectively recovered without any significant amount of 
defaults and vacancies. For the two solid layers just be- 
low the interface the situation is more complex. To dis- 
tinguish between different symmetries, we first perform 
a Voronoi construction for all the atoms in the layer. We 
then collect the set of first neighbors for each atom. 

For a fee solid with lattice parameter a, on the square 
lattice of the (100) orientation an atom has four nearest 
neighbors at a distance a/V2 and four second nearest 
neighbors at a distance a. On the other hand, for a tri- 
angular lattice (as the one of the (111) plane) the six 
neighbors all lie at the same distance a/y/2. In Fig. 6, 
we show a snapshot of the interface solid layer where the 
Delaunay triangulation is only drawn for the atoms that 
have six first neighbors at comparable distances, in or- 
der to reveal the local triangular structure. It is clear 
that most of the atoms have reached their positions on 
the square lattice but several islands with a triangular 
symmetry remain. Note that the number of atoms in 
the layer has already attained the value it will have deep 
in the solid with a perfect square structure. To compen- 
sate for the higher density of the triangular structure, the 
corresponding islands are surrounded by a border region 
where the density is very low. This coexistence of two 
symmetries is not observed in our smallest system: one 
can imagine that for a small area the square structure 
is easily formed and hence triangular islands do not ap- 
pear. This phenomenon is very close to the well known 
reconstruction of the (100) solid-vapor interface where 
the first layer adopts a triangular structure pl[ . Turning 
back to the solid-liquid interface, the system apparently 
uses some of the solidification driving force to eliminate 
one of the two phases and finally reach an almost perfect 
square symmetry. Hence, the interface velocity is lower 
for larger systems. 

Such an in-plane ordering is not taken into account 
in the collision-limited model but in spite of this we re- 
cover the predicted y/2 value for the ratio Hwo/Huo- This 
suggests a similar effect, roughly of the same order, for 
the (110) orientation. We have not been able to visual- 
ize ordering at (110) interfaces but one could imagine a 
mechanism reminiscent of the missing row reconstruction 
observed for (110) solid-vapor interfaces. 
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FIG. 6. Snapshot showing the atoms in the (100) solid 
layer next to the interface. The Delaunay triangulation is 
only drawn in regions with triangular underlying symmetry. 




FIG. 7. Velocity of the (111) solid-liquid interface as a 
function of undercooling. 



orientation. The extrapolated value of the kinetic coeffi- 
cient. 



Hm = 7.0±1.0cms~ 1 K- 1 , 



(12) 



is now 60 percent below its value for S — So. Relatively 
to the two other orientations, we find 



Hm ~ 0.37/iioo — 0.56/ino- 



(13) 



IV. THE SPECIAL CASE OF (111) INTERFACE 

We now turn to the case of the (111) orientation. In 
the same way as above for the (100) and (110) orienta- 
tions we calculate the interface temperature for different 
velocities. As can be seen in Fig. 7, a linear kinetic 
law is also valid for the (111) orientation. Results of the 
finite-size analysis, presented in Fig. 5, show that the 
size effects are much more pronounced than for the (100) 



These ratios largely differ from the values predicted by 
Eq. (1), respectively 2/^3 ~ 1.15 and 2t/2/3 ~ 1.63. 
The (111) orientation, expected to grow faster because 
of a larger interlayer spacing, is surprisingly found to be 
the slowest one. This discrepency tells us that the growth 
mechanism for the (111) orientation is not, or at least not 
only, a collision-limited one. 

Here again we look at the symmetries inside the lay- 
ers close to the interface. For a (111) layer, there are 
three possible ordered phases lying on three different but 
equivalent sub-lattices that we will call a, b and c. As the 
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FIG. 8. Snapshots of the three solid layers immediately 
below the inteface (top layer in contact with the liquid phase) . 
The grey levels correspond to the three different sub-lattices: 
a (white), b (light grey), c (dark grey). 



stacking fault energy is weak for gold (it is actually zero 
for the potential we use), once a perfect, say o, layer is 
formed, the next layer to form is either b or c. In Fig. 8, 
we show a snapshot of the three uppermost solid layers 
and we distinguish between atoms belonging to a, b or 
c phases. For the lowest solid layer, phase a is selected 
and it occupies the whole plane. For the layer just above, 
there is coexistence between b and c sub-lattices. Finally, 
in the highest solid layer, all three phases coexist. We re- 
cover here the effect first observed by Broughton et al. 
fllf for a LJ potential. For a (111) orientation the system 
hesitates between the different phases it can equivalently 
form. Here again, the system dissipates a part of the 
available driving force to select one of the phases. As a 
consequence, the velocity of the interface is reduced as 
compared to the value expected for a purely collision- 
limited growth. The size effect is easily understood be- 
cause in a small system, coexistence is strongly reduced. 
It would be of interest to determine the amount of driv- 
ing force spent in this in-plane organisation in order to 
estimate the corresponding decrease in Km. To perform 
this, one could for instance use a 3-state Potts model 
in three dimensions with ferromagnetic intra-plane and 
anti-ferromagnetic inter-plane interactions. To conclude 
this section we have to point out that phase coexistence 
is related to the value of the stacking fault energy E s . 
For a material with large E s phase coexistence should be 
less probable and the front velocity in better agreement 
with the prediction of Eq. (1). 



V. ASYMMETRY BETWEEN MELTING AND 
SOLIDIFICATION 

As discussed in the introduction, asymmetry is obvious 
for faceted materials but is not as clear when consider- 
ing rough materials like metals. The question is to know 
if, at equal absolute undercooling and superheating, the 
solid-liquid and liquid-solid fronts have the same velocity. 



With our simulation scheme, this study is straightfor- 
ward, since both a melting and a solidification fronts are 
simulated at once: no additional calculations are thus 
required. Fig. 9 represents the velocities of both the 
melting and solidification fronts as functions of T^ — Ti for 
the (111) orientation (in our conventions a positive veloc- 
ity corresponds to solidifation) . The data are obtained 
in a system of size S = Sq and corrected according to 
the finite-size analysis reported above. It is important to 
note that no size effects arc actually found for the melting 
front: in contrast with the solidification front, the melt- 
ing interface temperature remains the same whatever the 
system size. This can be understood if one remembers 
that for solidification, especially for the (111) orientation, 
growth is not only collision-limited but also requires in- 
plane ordering. This is no longer the case for melting, 
which justifies the absence of size effects. The asymme- 
try shown in Fig. 9 is larger for the (111) orientation. 
The same analysis is also made for the two other orien- 
tations and we find the following degrees of asymmetry: 
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where the superscripts s and m refer respectivly to solid- 
ification and melting kinetics. An asymmetry is revealed 
in the three cases but it is more pronounced for the (111) 
orientation in the same way as size effects observed dur- 
ing solidification. We conclude here that this asymmetry 
is directly related to the ordering within the interface 
layers. The asymmetry is strong for (111) because of 
the peculiar growth mechanism discussed in the previous 
section. 

The melting front is found to be faster than the solid- 
ification interface in agreement with the idea that disor- 
dering is an easier task than ordering. Our results con- 
firm the majority of experimental and numerical studies 
|]l4| , p2"| -|2"5f| . We also confirm the conclusions of a debate 
between Richards ||q] and Oxtoby and co-workers [^7|j28| 
on the importance of density change on the asymmetry 
between melting and solidification kinetics. In agreement 
with the conclusions of Oxtoby, the gold density change 
at melting is small (~ 2%) and can not be responsible 
for such an important asymmetry. On another hand, 
Tepper |l5f| does not find asymmetry for the growth of 
a (100) LJ solid. Even if the materials differ, they both 
belong to the same class of rough materials and such a 
qualitative difference may be surprising. Nevertheless, 
one should remember the strong tendency to surface re- 
construction in Au, as observed for the (100) orientation 
where triangular-like regions are formed. This tendency 
is futhermore enhanced by the use of Ercolessi glue poten- 
tial but is weaker for a LJ potential, what could explain 
the different behaviors observed. 
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FIG. 9. Velocity of the (111) solidification and melting 
interfaces as a function of undercooling. Note that the results 
for solidification incorporate finite-size corrections. 

Finally, comparing the melting kinetic coefficients in 
the different orientations, we find /Xj^ > /i^u > nTlo- 
We presently do not have a satisfactory explanation for 
this hierarchy in the melting kinetics. 

VI. DISCUSSION 

Our molecular dynamics simulations of zone melting 
experiments allow us to measure simultaneously the so- 
lidification and melting kinetics for a pure element. 

For (100) and (110) orientations, growth is apparently 
well described by a collision-limited process. Neverthe- 
less, we observe small 2D islands with triangular sym- 
metry to form in the solid layer at the (100) solid-liquid 
interface. As a consequence, size effects and asymmetry 
between melting and solidification are found. We can not 
decide whether this effect is solely due to the tendency 
of the glue potential to overestimate surface reconstruc- 
tion, or if it is an intrinsic property of gold and/or other 
metals. 

The case of the (111) orientation is rather special. 
Phase coexistence of three triangular sub-lattices, as first 
proposed by Broughton et al. |ll| , is recovered. This pe- 
culiar behavior has a strong influence on the kinetics of 
the interface. Our finite-size analysis show that in order 
to measure a realistic value of the kinetic coefficient one 
has to simulate systems with a solid-liquid interface area 
larger than 100 x 100A 2 . The consequence on asymmetry 
between melting and solidification is also of importance: 
for a given driving force, the melting front is more than 
three times faster than the solidification one. Because of 
this disagreement with a purely collision-limited growth, 
no analytical model seems, at present, able to predict 
the kinetic law of a (111) interface. As discussed previ- 
ously, it would be interesting to use a statistical model to 



extract the amount of driving force spent for phase sepa- 
ration in order to modify Eq. (1) and find an acceptable 
expression for the interface velocity. 

For melting we find the following order between the 
different kinetic coefficients : n^oo > n™ n > /x™ . To our 
knoweldge this hierarchy does not obey any existing law. 
This result will hopefully stimulate further investigations 
to reach a clear understanding of the specifities of melt 
growth as compared to crystal growth. 

The present study is devoted to the linear relationship 
between velocity and undercooling. For all the orienta- 
tions considered here, nonlinear effects appear at velocity 
V ~ 30 ms^ 1 and undercooling AT ~ 200i^. It is not 
possible to explain this deviation using either the diffu- 
sion limited |^,|3(J or the collision-limited growth law. 
This suggests a possible change in the interface structure 
for such large deviations from equilibrium. Density differ- 
ence between the liquid and solid phases should also con- 
tribute to trigger nonlinear behavior |26| . Understanding 
this cross-over would be of importance in the context of 
very rapid solidification. 

Finally, we would like to stress that the kinetic effects 
can contribute to the anisotropy of the segregation coef- 
ficient k(V) for a binary alloy. At sufficiently large veloc- 
ity, one expects an important difference in the interface 
temperatures for (111) and (100) orientations. As a con- 
sequence, the diffusivity of solvent atoms and hence the 
segregation coefficient, as predicted by the Aziz law [ MJ , 
should also differ. This effect may cause solute trapping 
to appear at lower velocities for (111) than for (100) or 
(110) orientations. We are currently investigating such 
segregation effects induced by kinetic anisotropy in the 
Al-Cu system. 
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